!================================================================
!                      v0(j),var(j)
! ------------#----------   ----- j ----------------#------------
!      fm(j-1)|fp(j-1)               fm(j)|fp(j)
!           flux(j-1)                 flux(j)
!
! fp(j): ---(-3)---(-2)---(-1)---(0)---(1)---(2)---(3)---(4)---
!         1: |--------------------|
!                 2: |------------------|
!                       3: |------------------|
!                              4: |------------------|
!
! fm(j): ---(-3)---(-2)---(-1)---(0)---(1)---(2)---(3)---(4)---
!                                    1: |-----------------|
!                              2: |-----------------|
!                       3: |------------------|
!                 4: |------------------|
!================================================================


subroutine wenozn7balsara_p( fps, fpx )
    use commondata
    implicit none
    
    real*8 :: fps(-3:3), fpx
    real*8 :: omega(4), fu(4), deriv(3,4), beta(4)
    real*8 :: tao7, tao12, ac, aa
    real*8 :: totalpha, alpha(4), fomega(4)

    omega(1) =  1.d0/35.d0;  omega(2) = 12.d0/35.d0
    omega(3) = 18.d0/35.d0;  omega(4) =  4.d0/35.d0
    
    fu(1) = -1.d0/4.d0*fps(-3) + 13.d0/12.d0*fps(-2) - 23.d0/12.d0*fps(-1) + 25.d0/12.d0*fps(0)
    fu(2) = 1.d0/12.d0*fps(-2) - 5.d0/12.d0*fps(-1) + 13.d0/12.d0*fps(0) + 1.d0/4.d0*fps(1)
    fu(3) = -1.d0/12.d0*fps(-1) + 7.d0/12.d0*fps(0) + 7.d0/12.d0*fps(1) - 1.d0/12.d0*fps(2)
    fu(4) = 1.d0/4.d0*fps(0) + 13.d0/12.d0*fps(1) - 5.d0/12.d0*fps(2) + 1.d0/12.d0*fps(3)
    
    deriv(1,1) = 1.d0/60.d0*( -19.d0*fps(-3) + 87.d0*fps(-2) - 177.d0*fps(-1) + 109.d0*fps(0) )
    deriv(2,1) = 1.d0/2.d0 *( -fps(-3) + 4.d0*fps(-2) - 5.d0*fps(-1) + 2.d0*fps(0) )
    deriv(3,1) = 1.d0/6.d0 *( -fps(-3) + 3.d0*fps(-2) - 3.d0*fps(-1) + fps(0) )
    deriv(1,2) = 1.d0/60.d0*( 11.d0*fps(-2) - 63.d0*fps(-1) + 33.d0*fps(0) + 19.d0*fps(1) )
    deriv(2,2) = 1.d0/2.d0 *( fps(-1) - 2.d0*fps(0) + fps(1) )
    deriv(3,2) = 1.d0/6.d0 *( -fps(-2) + 3.d0*fps(-1) - 3.d0*fps(0) + fps(1) )
    deriv(1,3) = 1.d0/60.d0*( -19.d0*fps(-1) - 33.d0*fps(0) + 63.d0*fps(1) - 11.d0*fps(2) )
    deriv(2,3) = 1.d0/2.d0 *( fps(-1) - 2.d0*fps(0) + fps(1) )
    deriv(3,3) = 1.d0/6.d0 *( -fps(-1) + 3.d0*fps(0) - 3.d0*fps(1) + fps(2) )
    deriv(1,4) = 1.d0/60.d0*( -109.d0*fps(0) + 177.d0*fps(1) - 87.d0*fps(2) + 19.d0*fps(3) )
    deriv(2,4) = 1.d0/2.d0 *( 2.d0*fps(0) - 5.d0*fps(1) + 4.d0*fps(2) - fps(3) )
    deriv(3,4) = 1.d0/6.d0 *( -fps(0) + 3.d0*fps(1) - 3.d0*fps(2) + fps(3) )
    
    beta(:) = ( deriv(1,:)+0.1d0*deriv(3,:) )**2 + 13.d0/3.d0*deriv(2,:)**2 + 781.d0/20.d0*deriv(3,:)**2
    
    tao7 = abs( beta(1) - beta(4) )
    tao12 = ( fps(-3) - 6.d0*fps(-2) + 15.d0*fps(-1) - 20.d0*fps(0) &
            + 15.d0*fps(1) - 6.d0*fps(2) + fps(3) )**2
    !aa = 10.d0
    ac = ( ( beta(1) + beta(4) - tao7 + eps )/( tao7 + eps ) )**2
    
    alpha(:) = omega(:)*( ac + ( tao12/(beta(:)+eps) )**2 )
    !alpha(:) = omega(:)*( ac + ( (tao12+10.d0*tao7**2)/(beta(:)+eps) )**2 )
    totalpha = alpha(1) + alpha(2) + alpha(3)+ alpha(4)
    fomega(:) = alpha(:)/totalpha
    fpx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)+ fu(4)*fomega(4)
    
    return
    end subroutine wenozn7balsara_p


subroutine wenozn7balsara_m( fms, fmx )
    use commondata
    implicit none
    
    real*8 :: fms(-2:4), fmx
    real*8 :: omega(4), fu(4), deriv(3,4), beta(4)
    real*8 :: tao7, tao12, ac, aa
    real*8 :: totalpha, alpha(4), fomega(4)

    omega(1) = 1.d0/35.d0;   omega(2) = 12.d0/35.d0
    omega(3) = 18.d0/35.d0;  omega(4) = 4.d0/35.d0
    
    fu(1) = -1.d0/4.d0*fms(4) + 13.d0/12.d0*fms(3) - 23.d0/12.d0*fms(2) + 25.d0/12.d0*fms(1)
    fu(2) = 1.d0/12.d0*fms(3) - 5.d0/12.d0*fms(2) + 13.d0/12.d0*fms(1) + 1.d0/4.d0*fms(0)
    fu(3) = -1.d0/12.d0*fms(2) + 7.d0/12.d0*fms(1) + 7.d0/12.d0*fms(0) - 1.d0/12.d0*fms(-1)
    fu(4) = 1.d0/4.d0*fms(1) + 13.d0/12.d0*fms(0) - 5.d0/12.d0*fms(-1) + 1.d0/12.d0*fms(-2)
    
    deriv(1,1) = 1.d0/60.d0*( -19.d0*fms(4) + 87.d0*fms(3) - 177.d0*fms(2) + 109.d0*fms(1) )
    deriv(2,1) = 1.d0/2.d0 *( -fms(4) + 4.d0*fms(3) - 5.d0*fms(2) + 2.d0*fms(1) )
    deriv(3,1) = 1.d0/6.d0 *( -fms(4) + 3.d0*fms(3) - 3.d0*fms(2) + fms(1) )
    deriv(1,2) = 1.d0/60.d0*( 11.d0*fms(3) - 63.d0*fms(2) + 33.d0*fms(1) + 19.d0*fms(0) )
    deriv(2,2) = 1.d0/2.d0 *( fms(2) - 2.d0*fms(1) + fms(0) )
    deriv(3,2) = 1.d0/6.d0 *( -fms(3) + 3.d0*fms(2) - 3.d0*fms(1) + fms(0) )
    deriv(1,3) = 1.d0/60.d0*( -19.d0*fms(2) - 33.d0*fms(1) + 63.d0*fms(0) - 11.d0*fms(-1) )
    deriv(2,3) = 1.d0/2.d0 *( fms(2) - 2.d0*fms(1) + fms(0) )
    deriv(3,3) = 1.d0/6.d0 *( -fms(2) + 3.d0*fms(1) - 3.d0*fms(0) + fms(-1) )
    deriv(1,4) = 1.d0/60.d0*( -109.d0*fms(1) + 177.d0*fms(0) - 87.d0*fms(-1) + 19.d0*fms(-2) )
    deriv(2,4) = 1.d0/2.d0 *( 2.d0*fms(1) - 5.d0*fms(0) + 4.d0*fms(-1) - fms(-2) )
    deriv(3,4) = 1.d0/6.d0 *( -fms(1) + 3.d0*fms(0) - 3.d0*fms(-1) + fms(-2) )

    beta(:) = ( deriv(1,:)+0.1d0*deriv(3,:) )**2 + 13.d0/3.d0*deriv(2,:)**2 + 781.d0/20.d0*deriv(3,:)**2
    
    tao7 = abs( beta(1) - beta(4) )
    tao12 = ( fms(-2) - 6.d0*fms(-1) + 15.d0*fms(0) - 20.d0*fms(1) &
            + 15.d0*fms(2) - 6.d0*fms(3) + fms(4) )**2
    !aa = 10.d0
    ac = ( ( beta(1) + beta(4) - tao7 + eps )/( tao7 + eps ) )**2
    
    alpha(:) = omega(:)*( ac + ( tao12/(beta(:)+eps) )**2 )
    totalpha = alpha(1) + alpha(2) + alpha(3)+ alpha(4)
    fomega(:) = alpha(:)/totalpha
    fmx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)+ fu(4)*fomega(4)
    
    return
    end subroutine wenozn7balsara_m